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ABSTRACT 

We consider the information that can be derived about massive black-hole binary populations 
and their formation history solely from current and possible future pulsar timing array (PTA) 
results. We use models of the stochastic gravitational-wave background from circular massive 
black hole binaries with chirp mass in the range 10^ — 10* W© evolving solely due to radiation 
reaction. Our parameterised models for the black hole merger history make only weak assump¬ 
tions about the properties of the black holes merging over cosmic time. We show that current 
PTA results place an upper limit on the black hole merger density which does not depend on 
the choice of a particular merger history model, however they provide no information about 
the redshift or mass distribution. We show that even in the case of a detection resulting from a 
factor of 10 increase in amplitude sensitivity, PTAs will only put weak constraints on the source 
merger density as a function of mass, and will not provide any additional information on the 
redshift distribution. Without additional assumptions or information from other observations, a 
detection cannot meaningfully bound the massive black hole merger rate above zero for any 
particular mass. 
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1 INTRODUCTION 

Massive black holes (MBHs) reside at the centre of most galaxies 
(see e.g. Kormendy & Ho 2013, and references therein), and are 
believed to have a central role in their evolution (see e.g. Volonteri 
2012, and references therein for a recent review). Mapping the 
population of MBHs, studying their properties, demographics, and 
their connection to the broader formation of structure is one of the 
open problems of modern astrophysics. This is however difficult 
to tackle, due to the large range of scales and the wide variety of 
physical processes involved. The MBH evolutionary path remains 
a highly debated subject with many competing hypotheses still in 
play. Currently favoured hierarchical structure formation scenarios 
imply frequent galaxy mergers (White & Rees 1978). As a result 
MBH binaries (MBHBs) should be quite common in the Universe 
(Begelman et al. 1980; Volonteri et al. 2003). To-date there is no 
confirmed observed MBHB, although a number of candidates exist 
(see e.g. Dotti et al. 2012, and references therein) and tantalising 
claims have been recently made (Graham et al. 2015b; Liu et al. 
2015; Graham et al. 2015a). 

A means to survey MBHBs is through the observation of gravi¬ 
tational waves (GWs) that these systems generate as they inspiral 
towards their final merger. The accurate timing of an array of highly- 
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stable millisecond pulsars - a Pulsar Timing Array (PTA, Foster 
& Backer 1990) - provides a direct observational means to probe 
the cosmic population of MBHBs on orbital timescales of order 
of several years. Astrophysical modelling suggests that the radia¬ 
tion emitted by an ensemble of MBHBs produces a GW stochastic 
background in the frequency range ~ 10^® — 10^^ Hz, where PTAs 
operate (Sesana et al. 2008, 2009; Ravi et al. 2012; Sesana 2013b). 
Such a background affects the time of arrival of radio pulses in a 
characteristic fashion (Sazhin 1978; Detweiler 1979; Hellings & 
Downs 1983), which can be used to discriminate the signal from a 
plethora of other undesired effects (Lentati et al. 2015). 

Over the last decade pulsar timing has been used to put progres¬ 
sively tighter constraints on gravitational radiation in this frequency 
regime, (see, e.g. Jenet et al. 2006). More recently the three in¬ 
ternational consortia consisting of the Parkes Pulsar Timing Array, 
PPTA (Shannon et al. 2015; Shannon et al. 2013), NANOGrav (De- 
morest et al. 2013; Arzoumanian et al. 2015) and the European 
Pulsar Timing Array, EPTA (Lentati et al. 2015; van Haasteren et al. 
2012, 2011), which in collaboration form the International Pulsar 
Timing Array (Hobbs et al. 2010), have used data from observations 
of unprecedented sensitivity to place constraints that are starting 
to probe astrophysically interesting regions of the parameter space 
(Sesana 2013b). 

In this paper, we consider a GW stochastic background pro¬ 
duced by MBHBs in circular orbits losing energy and angular mo- 



L2 H. Middleton et al. 


mentum purely through GW emission. We use an analytical merger 
rate model which makes minimal assumptions about the cosmologi¬ 
cal history of MBHB evolution and can capture the key characteris¬ 
tics of simulation results to investigate the astrophysical implications 
of current (Shannon et al. 2013; Lentati et al. 2015; Arzoumanian 
et al. 2015), and future plausible (Siemens et al. 2013; Ravi et al. 
2015) PTA results (either an upper-limit or a detection). Because our 
model is fully general—not committing to any particular cosmolog¬ 
ical MBHB merger history—we can identify and separate features 
of the merger history that are constrained by PTA data alone, from 
those that can only be constrained by adopting a particular merger 
history (as e.g. done in Shannon et al. (2013) and Arzoumanian 
et al. (2015))—in other words, by applying a particular cosmologi¬ 
cal prior. Because our model is capable of reproducing the MBHB 
cosmic population found in cosmological simulations for certain 
choices of parameters, our results will be consistent with (but much 
broader than) those that would be obtained under a choice of specific 
classes of MBHB merger history models. 

In Section 2 we describe our method and the model used for 
the merger rate. In Section 3 we present our results for several 
upper limits and for a possible future detection, and we discuss 
their implications for the population of MBHBs. We present our 
conclusions in Section 4. 


2 MODEL AND METHOD 
2.1 Astrophysical model 

For the standard scenario of circular binaries driven by radiation 
reaction only, the characteristic strain of the GW stochastic back¬ 
ground, h at frequency / is (Phinney 2001): 

= 3^^/3g2 j ^ 

d^N 

dVYzdlogiQ^ ’ 

where z is the redshift and is the chirp mass related to the binary 
component masses (mi, m 2 ) by = (mim 2 )^/^j(m\ +m 2 )*/^. 
The integral sums over the sources in z and weighted by the 
distribution of the source population, YN/dVYzdlogiQ,y£, the 
number of binary mergers per co-moving volume, redshift and (rest- 
frame) chirp mass interval. We choose a simple model for this, 
described by 

d'^N _ ■ 

dVcdzdlogiQj^ IO’Mq 

X [(l+z)'^exp-(^/®)] (2) 

where is the time in the source rest-frame (here we use Hq = 
70km s^*Mpc^', SIm = 0.3, IIa = 0-7 and = 0). Following 
general astrophysical assumptions, we consider a scenario where the 
GW background is produced by MBHBs in the redshift and chirp 
mass range of 0 < z < 5 and 10^ < „#/M© < 10**. These ranges 
set the integration limits of Eqn. (1). 

The model is described by five parameters. The parameter hq 
is the normalised merger rate per unit rest-frame time, co-moving 
volume and logarithmic (rest-frame) chirp mass interval. The pa¬ 
rameters j8 and zo describe the distribution of the sources in redshift. 
The parameter j3 controls the low-redshift power-law slope and the 
parameter zo the high-redshift cut-off for the distribution; the peak 
of the merger rate d^N/dt^dV^, corresponds to a redshift (zoj8 — 1). 



The parameters a and .y#* provide a similar description of the 
chirp mass distribution. The model was chosen to capture the ex¬ 
pected qualitative features of the cosmic MBH merger rate without 
restricting to any particular merger history; for example, it can re¬ 
produce rates extracted from merger tree models (Volonteri et al. 
2003; Sesana et al. 2008), and large scale cosmological simulations 
of structure formation (Springel et al. 2005; Sesana et al. 2009). 

The characteristic amplitude has a simple power-law scaling, 
and we can re-write Eqn. (1) as 

where Aiyf is the characteristic amplitude at the reference frequency 
/lyr = lyr^*, which is customarily used when quoting limits in 
the PTA literature. A single number, the amplitude Ajy^, carries the 
whole information about the merging history of MBHBs (within the 
model considered in this paper), that one wishes to reconstruct from 
the observations. 


2.2 Method 


The objective is to put constraints on the population parameters, 
which we denote by 6, given the results of PTA analyses. In our 
case 0 is a 5-dimensional parameter space, 6 = {ho,l5,zo, a,..#*}. 
We want to compute the posterior density function (PDF) of the 
parameters given PTA observations denoted by d. The population pa¬ 
rameters fully specify the gravitational wave signal /;(/; 0) (Eqs. (1) 
and (3)), which in turn specifies the statistical properties of the GW- 
induced deviations to pulse arrival times, the PTA observable. Given 
data from pulsar timing and our model for the merger rate (Eqn (2)), 
we use Bayes’ theorem to find the posterior distribution of the model 
parameters p{6\d). 




(4) 


where p(d|Aiyr(6)) is the PTA likelihood for a given stochastic 
background, h(f\Q),p{8) is the prior on the model parameters and 
p(d) is the evidence. In standard analysis of the PTA data, con¬ 
straints are put on the GW characteristic amplitude at periods of 
one year, Aiyr, which in turn is a function of the parameters of 
the underlying population, specified by h{f\ 0). The PTA analysis 
uses a likelihood function, p{d\Aiyy{6)), which we approximate 
as described below. Our method does not rely on this approxima¬ 
tion; we use it only for analytical convenience in this paper. If a 
given PTA analysis provides a posterior distribution for Ajy^ then a 
straightforward re-weighting can produce the corresponding likeli¬ 
hood required for our analysis (if flat priors on Aiyj are used in the 
analysis then the re-weighting is trivial because the posterior and 
the likelihood are proportional to each other). 

In this paper we consider the two cases in which the PTA 
analysis provides either an upper-limit or a detection. For the upper- 
limit scenario we model p{d\Aiyd using a Fermi-like distribution: 


Pul(rf|Alyr)°= (exp((Aiyr-Aui)/(Tui) + l) (5) 


where Aui is the upper-limit value returned by the actual analysis 
and the sharpness of the tail-off, CTui can be adjusted to give an upper 
limit with a chosen confidence, which we set at 95%. We model a 
detection scenario using a Gaussian in the logarithm of Ajy^: 


Pdet (rf|Alyr) exp(-(logio(Alyr) - logio(Adet))V20’|et) (6) 
at a chosen level of detection, Ajet- We choose the width of the 
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Figure 1. Posteriors for the merger rate density. The top row shows the merger rate density in chirp mass (integrated over redshift), (f-N/dVcdiog^Q^ and the 
bottom row in redshift (integrated over chirps mass), d^N/dVcdz, for two 95% confidence upper limits at 1 x 10“^^ (left) and 1 x 10“^^ (centre) and a detection 
at 1 X 10“^^ (right), as described in the text. We consider contributions to the gravitational wave background from massive black holes in the chirp mass range 
10^ < <10^^ and redshift range 0 < z < 5. The solid black line gives the posterior median; dark grey, mid-grey and light-grey bands show the central 

68%, 95%, and 99% credible interval, respectively. The dashed lines show draws from the posterior. For comparison, the overlaid dark areas represent the 99.7% 
confidence regions predicted by the MBH assembly models of Sesana (2013b). For these models, we show only the chirp mass in the range w 10^Mq — 10®Mq 
as outside this interval the lower percentile is zero. For the redshift range, these models only consider MBHB mergers for z < 1.3. 


detection to be Odet = 0'2- We compute the marginalised distribu¬ 
tion on the model parameters 6 using two independent sampling 
techniques, to verify the results of our analysis: a nested sampling 
approach (Veitch & Vecchio 2010) and emcee, an ensemble Markov 
Chain Monte Carlo sampler (Foreman-Mackey et al. 2013). 

Our priors on the model parameters are set as follows. We use 
a prior on hp that is flat in logjghp down to a lower limit, which 
we set to Mo = 10 ^^®Mpc^^Gyr^*, after which it is flat in no to 
zero. The prior upper-bound is set to lO^Mpc^^Gyr^'. This value 
is set by the ultra-conservative assumption that all the matter in the 
Universe is formed by MB FIs. Our prior allows for the number of 
mergers to span many orders of magnitude (flat in log) but avoids 
divergence as hg —> 0. It also allows for the absence of MBHB 
binaries merging within an Hubble time. The priors for the other 
parameters are uniform within ranges that incorporate values that 
give a good fit to semi-analytical merger tree models (see e.g. Sesana 
etal. 2008, 2009; Sesana 2013b): a G [-3.0,3.0], e [-2.0,7.0], 
zo G [0.2,5.0] andlog[o,-#*/M 0 G [6.0,9.0]. While our prior allows 
for parameter values that can reproduce the merger rates of detailed 
models, it is uninformative in that we do not assume that the merger 
rate distribution must take values from those models. Our priors 
reflect large theoretical uncertainties about MBHB formation and 
evolution scenarios, and the lack of any confirmed MBHB candi¬ 
date, (see however Graham et al. 2015b; Liu et al. 2015; Graham 
et al. 2015a). 


Our method is summarised as: (i) produce a likelihood for 
Aiyr (in the case of an actual analysis by using smoothed posterior 
samples from PTA results, re-weighted if necessary depending on 
the prior), (ii) choose a model for the merger rate of MBHBs, (iii) 
produce posterior density functions for the model parameters from 
which we can infer properties of the MBHB population. 


3 RESULTS 

Current upper-limits on the GW stochastic background obtained 
recently are = 1 x 10^*^, 1.5 x 10^^^, 3 x 10^^^ for the 

PPTA (Shannon et al. 2015), NANOGrav (Arzoumanian et al. 2015), 
and the EPTA (Lentati et al. 2015) respectively. The sensitivity gain 
provided by the addition of new pulsars to the PTAs and more recent 
data sets may allow in the short-to-mid term to reach a sensitivity 
below Aiyr = 1.0 X 10^*^, and in the more distant future Aiyr ~ 
10^*® (Siemens et al. 2013; Ravi et al. 2015). As a consequence, 
here we consider three PTA analysis outcomes; (i) an upper-limit at 
95% confidence of 1 x 10^^^, which represents the present state of 
play and either (ii) an upper-limit (at 95% confidence) of 1 x 10^*® 
or (iii) a detection at the same level, that is A^et = 1 x 10 ^^®, and 
(7det = 0-2 in Eqn. ( 6 ), which describes possible results coming 
from the expected improvements of the PTA sensitivity in the next 
five-to-ten years. 
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The main results of our analysis are summarised in Fig. 1, 
which shows the inferred posterior distribution of the merger his¬ 
tory of MBHBs in terms of the MBHBs co-moving volume merger 
density per redshift and (logarithm of) chirp mass intervals. Fig. 2 
provides PDFs on selected parameters based on current PTA lim¬ 
its, and Fig. 3 provides a similar summary for a future limit or a 
detection at the level described above. 

We consider first the implications of current limits. The PDFs 
on the parameters mq and •^* of the model are shown in Fig. 2; 
we do not provide the equivalent plots for a, j3 and zo as they 
are equivalent to the prior. Fig. 2 clearly shows that the present 
PTA limits enable us to reduce the allowed normalization of the 
MB FIB merger rate density to ho iS 5 x lO^^Mpc^^Gyr^* with 
95% confidence, but yield no additional constraints on the other 
parameters of the model. 

Our model contains parameters describing the shape of the 
merger rate distribution in redshift and chirp mass. The PDFs of 
those parameters induce a posterior density on d^N/dVcdlogiQ ^ 
and d^N/dVcdz, integrating over redshift and chirp mass respec¬ 
tively, shown in Fig. 1. We see that current observations limit the 
maximum merger density as a function of mass, but place no con¬ 
straints on the shape of the distribution. The corresponding num¬ 
ber of sources per frequency bin that contribute to the signal is 
dN/dfdlogiQ^ and we find that for masses above a 

few X IO^Mq, our upper limit on d^N/dVcdlog^Q^ implies that 
at a frequency around 1.8 nPIz there is fewer than one source per 
frequency bin (taken to be A/ = 1 /T, with T = 17.66 yr, the times- 
pan of current EPTA datasets (Lentati et al. 2015)). This means 
that at those large masses, the assumption that the observed GW 
signal is stochastic is violated, and our analysis cannot be used to 
constrain the exact shape of the mass function here (in this case a 
different PTA search approach would be necessary, see e.g. Babak 
et al. (2015); Arzoumanian et al. (2014); Taylor et al. (2015)). While 
current PTA observations provide feeble constraints on the shape of 
the mass distribution, they yield no information about the redshift 
distribution. The bottom-left panel of Fig. 1 shows no structure in 
d^N/dVcdz. 

It is useful to compare these results to limits on the MBHB 
merger rates implied by binary candidates reported in the litera¬ 
ture and to specific theoretical models. Let us first consider what 
is known observationally today. A few MBHB candidates have 
been reported recently. Graham et al. (2015b) suggested the pos¬ 
sible observation of a MBHB at redshift z = 0.2784 with (rest 
frame) total mass \og{M/M q) ~ 8.5 and period of ~ 1884 days. 
Liu et al. (2015) reported the observation of a potential MBHB at 
z = 2.060 with a shorter period of 542 days and primary MBHB 
mass \og{M/M q) ~ 9.97. Using the redshift to calculate the en¬ 
closed volume and the binary parameters for the time to merger 
we can estimate the predicted rate from each of these observations. 
Assuming that these two systems are indeed MBHBs, and that their 
constituents are of comparable mass, they imply merger rates of 
3 X lO^^Mpc^^Gyr^* and k:; O.lMpc^^Gyr^^ (this latter num¬ 
ber takes into account that the source has been found in an analysis 
of only one of the Pan-STARRS 1 Medium Deep Survey fields of 
8 deg^). In turn they yield a merger density d^N/dV^dlogiQ^ ~ 
lO^^Mpc^^ at 3 X 10 *Mq and IMpc^* ^ 10*® Mq, 
respectively. The upper left panel of Figure 1 clearly shows that 
the rate density inferred from Graham et al. (2015b) is consistent 
with current upper limits, while that inferred from Liu et al. (2015) 
is several orders of magnitude above the 99% credible interval 
implied by current PTA results. It is therefore unlikely that this 
source is a MBHB with the claimed parameters. Other proposed 



Figure 2. Marginalised posterior distributions for selected astrophysical 
parameters for the case of 95% upper-limit of 1 x 10“*^, which corresponds 
to the current status of the observations. The marginalised PDF on the merger 
rate parameter, ho is shown in the left panel, and the marginalised PDF on 
(.^,,ho) in the right panel, where the contours mark the 67% (solid) and 
95% (dashed) confidence regions. In the left panel, the dashed lines mark 
the 95% confidence width (—20.8 < log[ho/Mpc“^Gyr“'] < —2.3) while 
the dotted line marks the 95% upper limit (log[ho/Mpc“^Gyr“’] < —3.3). 
The left hand side of the distribution in log(ho/Mpc“^ Gyr“') follows our 
prior, while the right hand side is determined by the PTA upper limit. 


MBHBs in the literature (Valtonen et al. 2012; Kun et al. 2014) 
imply merger density estimates of 5 x 10^^ Mpc^* at ^ ~ 10®M© 
and 3 x 10^®Mpc^* zA ^ ^ 3.5 x 10 *Mq, which are consistent 
with the current upper limits. 

On the theoretical side, current limits are consistent with the 
assumption that most Milky-Way-like galaxies contain a MBH in 
the mass range considered here that undergoes ~ 1 major merger 
in an Hubble time. The density of Milky-Way-like galaxies is 
lO^^Mpc^*, which yields an estimate of d^NjdVcdXogiQ^ ~ 
10^*Mpc^*, which is consistent with our results at ^ ~ IQ^Mq, 
appropriate for a typical MBHB forming in the merger of 
Milky-Way-like galaxies. We also compare the limits on the 
d^N/dVcdlogiQ^ and d^N/dV^dz with specific distributions ob¬ 
tained from predictions of astrophysical models for the cosmic 
assembly of MBHs. We consider the models presented in Sesana 
(2013b), extended to include the most recent MBH-galaxy scaling 
relations (Kormendy & Ho 2013). These models produce a central 
99% interval of Aiyr e [2 x 10^*®,4 x 10^*^]. The 99.7% confi¬ 
dence region in the merger density from those models is marked 
by a dark-shaded area in each panel of Fig. 1. Two conclusions 
can be drawn: (i) present MBHB population models are consistent 
with current PTA limits; (ii) those models are drawn from a very 
restricted prior range of the parameters that control the evolution of 
MBHBs, driven by specific assumptions on their assembly history. 
For example, in those models there is a one-to-one correspondence 
between galaxy and MBH mergers. Our results are consistent with 
the conclusions drawn by Shannon et al. (2013) about the implica¬ 
tions of the PPTA limit for the MBHB merger history. However, 
since Shannon et al. (2013) consider specific models that lie close to 
the upper end of the 99% credible range allowed by current limits, 
they emphasise the fact that PTA limits might soon be in tension 
with those specific classes of models. 

We turn now to consider what we could infer about the MBHB 
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merger history in the future as PTA sensitivity increases. For def¬ 
initeness we consider both an upper limit and a detection at the 
level of Aiyr = 10^^®. Selected marginalised PDFs on the model 
parameters are shown in Fig. 3, where we see a slight correlation in 
the 2-D marginalised PDF of (..#*,mo)> as expected. This is simply 
explained by considering the Schechter-like mass prohle of Eqn. 2: 
as the characteristic mass decreases, and therefore the exponen¬ 
tial cut-off of MB FIB progressively depletes the high-mass portion 
of the population, a given value of the GW characteristic amplitude 
allows for a larger overall normalisation, «o- The posterior MBFIB 
merger densities per logarithm of chirp mass and redshift are shown 
in Fig. 1. For the case of an upper-limit the results are qualitatively 
similar to the case of the present PTA upper-limit, simply scaled ac¬ 
cordingly. In particular, despite the much tighter limit on the overall 
merger rate we are still unable to place any meaningful constraint 
on the redshift distribution of merging MBFIBs. The overall merger 
density as a function of redshift shifts by two orders of magnitude 
and the same is true for the merger density as a function of mass. 
Note that a non-detection at this level might pose a serious challenge 
to currently favoured theoretical MBH assembly models with simple 
black hole dynamics, as shown in the upper-centre panel of Fig. 1. 

In the case of detection the posterior on the shapes of the merger 
rate distribution in redshift and chirp mass are plotted on the right 
panels of Fig. 1. We still obtain essentially no bounds on the shape 
of the merger rate density in redshift. We also obtain no meaningful 
lower bound on the merger rate density for chirp masses. That is, 
there is no chirp mass at which we can bound the merger density 
above a rate physically indistinguishable from zero; we know that 
some MBFIBs merge, but we cannot determine which ones. Addi¬ 
tional information, such as theoretical assumptions, electromagnetic 
observations constraining the mass spectrum of merging black holes 
(like those discussed earlier in this Section), or gravitational wave 
observations that measure the binary mass spectrum directly (such 
as those of an eLISA-like instrument (Amaro-Seoane et al. 2013)), 
are required to place any constraints on the masses of the merging 
systems. For example, if we accept the priors provided by Sesana 
(2013b), the mass function of merging MBFIBs can be determined 
more precisely, as shown by the overlap between our posterior and 
the dark band in the upper-right panel of Fig. 1. 


4 CONCLUSIONS 

We have considered the implications of current PTA limits on the 
GW stochastic background to constrain the merger history of MB- 
HBs. Using a general model for the mass and redshift evolution of 
MBFIBs in circular orbit driven by radiation reaction, we find that 
existing PTA results alone place essentially no constraints on the 
merger history of MBFIBs. We also find that even with an increase 
in amplitude sensitivity of an order of magnitude, and assuming that 
a detection is made, no bounds can be put on the functional form 
of the merger rate density in redshift and chirp mass unless addi¬ 
tional information coming through a different set of observations is 
available. 

Finally we want to caution the reader that the results presented 
here apply only within the model assumptions that have been made. 
We have considered a generic (and well justified) functional form 
for the MBFIB merger rate density, but if one chooses a significantly 
different form (and associated priors for the parameters), results 
could be different (even radically). Moreover, it has been suggested 
that physical effects other than radiation reaction, such as gas and/or 
interactions with stars (e.g. Kocsis & Sesana 2011; Sesana 2013a; 



Figure 3. Posterior distribution for the upper limit (left) and detection (right) 
at 1 X 10“**. The top panels show the one dimensional posterior distri¬ 
bution for the merger rate parameter, «o. The dashed lines mark the 95% 
confidence width (upper limit: —20.9 < logio(ho/Mpc ^Gyr *) < -4.2; 
detection: —11.7 < log[Q(ho/Mpc“^ Gyr“*) < —1.3) and the dotted line the 
95% upper limit (upper limit: log[Q(fio/Mpc“^Gyr“*) = —5.0; detection: 
log[o(ho/Mpc“^ Gyr“*) = —1.9). The central and bottom panels show the 
two dimensional posterior distributions for ho with the mass parameters a 
and The solid and dashed contours mark the 67% and 95% confidence 
regions respectively. 


Sampson et al. 2015), could affect the evolution of MBHBs. These 
effects are not included in our model, and their impact on astrophys- 
ical inference needs to be evaluated in the future. 
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